RVAideMemoireパッケージにあるfisher.multcomp関数で,p.adjust関数にあるp値調整法でFisher検定の多重比較ができる.

データ例

library(RVAideMemoire)
d <- matrix(c(17,23,12,24,20,10),ncol=2,
            dimnames=list(c("cont","trt1","trt2"),
                          c("alive","dead")),
            byrow=TRUE)
d
##      alive dead
## cont    17   23
## trt1    12   24
## trt2    20   10

fisher.test関数

データ例で,conttrt1trt2の3群の生存率の比較をFisher検定で行う.

res12 <- fisher.test(d[c(1,2),])
res13 <- fisher.test(d[c(1,3),])
res23 <- fisher.test(d[c(2,3),])
c(`cont vs trt1`=res12$p.value,
  `cont vs trt2`=res13$p.value,
  `trt1 vs trt2`=res23$p.value)
## cont vs trt1 cont vs trt2 trt1 vs trt2 
##   0.48195094   0.05567911   0.01279792

RVAideMemoireパッケージのfisher.multcomp関数で, p.method="none"で調整なしのFisher検定のp値.

fisher.multcomp(d, p.method = "none")
## 
##         Pairwise comparisons using Fisher's exact test for count data
## 
## data:  d
## 
##         cont   trt1
## trt1 0.48195      -
## trt2 0.05568 0.0128
## 
## P value adjustment method: none

検定の過誤

\(m\)個の帰無仮説\(H_1,H_2,\ldots,H_m\)の検定をする状況を考える. \(m_0\)個の帰無仮説には真に正しい仮設とする. 下の表は,\(V\)個が第1種の過誤(真である帰無仮説を誤って棄却)が生じた検定数であり, \(R\)は棄却した検定数である. \(R\)は観測可能な確率変数に対し,\(S\), \(T\), \(U\), \(V\)は観測できない確率変数である. \(m\)\(m_0\)は定数であり,\(m_0\)は未知である.

\[\begin{array}{cccc} \hline {\rm Hypothesus} & {\rm Not\ Rejected} & {\rm Rejected} & {\rm Total} \\ \hline {\rm True} & U & V & m_0 \\ {\rm False} & T & S & m-m_0 \\ \hline {\rm Total} & m-R & R & m \\ \hline \end{array}\]

個々の帰無仮説\(H_i\)に対する検定における第1種の過誤は有意水準\(\alpha\)で制御されていが, \(m\)個の帰無仮説の検定の全体での第1種の過誤は,FWER(Family Wise Error Rate)で考える.

\[ {\rm FWER}=\Pr(V>0) \]

次に,棄却した検定数\(R\)に含まれる偽陽性(False Positive)の検定数として\(V\)を扱う方法を考える. 間違って棄却された帰無仮説の割合を\(Q=V/R\)とする.ただし,\(R=0\)のときは\(Q=0\)とする. この\(Q\)も確率変数である.その期待値をFDR(False Discovery Rate)と呼ぶ. \[ {\rm FDR}=E(Q)=E\left(\frac{V}{R} \right) \]

FWERとFDRには以下の関係がある.

ボンフェローニ法

\(m\)個の検定の\(p\)値を値が小さい順に並べ替えた\(p_1<p_2< \cdots <p_m\)に対して, ボンフェローニ法では,調整\(p\)値を\(q_i=\min(1,mp_i)\) \((i=1,2,\ldots,m)\)とする. 上の例だと,\(p_1=0.0128\)\(p_2=0.05568\)\(p_3=0.48195\)より,

\[q_1=3\times 0.0128=0.0384, \quad p_2=3\times 0.05568=0.1670, \quad p_3=3\times0.48195 \to 1\]

fisher.multcomp関数でp.method="bonferroni"とする.

fisher.multcomp(d, p.method = "bonferroni")
## 
##         Pairwise comparisons using Fisher's exact test for count data
## 
## data:  d
## 
##       cont    trt1
## trt1 1.000       -
## trt2 0.167 0.03839
## 
## P value adjustment method: bonferroni

ホルム法

ホルム法では次のように\(p\)値を調整し,\(q_i=(m-k+1)p_{k}<\alpha\)となる最大の\(k\)を見つけ,第1順位から第\(k\)順位までの帰無仮説を棄却し,第\((k+1)\)順位以降の帰無仮説の判定を保留する.

\[ q_1=mp_1,\ q_2=(m-1)p_2,\ q_3=(m-2)p_2,\ \cdots \ , q_{m-1}=2p_{m-1},\ q_m=p_1,\ \]

上の例だと,\(p_1=0.0128\)\(p_2=0.05568\)\(p_3=0.48195\)より,

\[q_1=3\times 0.0128=0.0384, \quad p_2=2\times 0.05568=0.11136, \quad p_3=1\times 0.48195=0.48195\]

fisher.multcomp関数でp.method="holm"とする.

fisher.multcomp(d, p.method = "holm")
## 
##         Pairwise comparisons using Fisher's exact test for count data
## 
## data:  d
## 
##        cont    trt1
## trt1 0.4820       -
## trt2 0.1114 0.03839
## 
## P value adjustment method: holm

FDR法

Benjamini-Hochberg法では次のように\(p\)値を調整し,\(q_i=(m/i)p_i < \alpha\)となる最大の\(k\)を見つけ, \(H_1,H_2,\ldots,H_k\)の仮設を棄却する. \(q_1<\alpha\)を満たさなかった場合は,どの帰無仮説も棄却しない.

\[ q_1=\frac{m}{1}p_1,\ q_2=\frac{m}{2}p_2,\ q_3=\frac{m}{3}p_2,\ \cdots \ , q_{m-1}=\frac{m}{m-1}p_{m-1},\ q_m=\frac{m}{m}p_m,\ \]

fisher.multcomp(d, p.method = "BH")
## 
##         Pairwise comparisons using Fisher's exact test for count data
## 
## data:  d
## 
##         cont    trt1
## trt1 0.48195       -
## trt2 0.08352 0.03839
## 
## P value adjustment method: BH